Fit condition model with environmental and biological predictors

Fit main model (see exploratory scripts and model comparison), visualize results.

Load packages

library(tidyverse); theme_set(theme_light(base_size = 12))
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(viridis)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(gganimate)
library(gifski)
library(latex2exp)
library(patchwork)
library(png)
library(RCurl)
library(wesanderson)
library(qwraps2) 
library(ggcorrplot)
library(ggridges)
# remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)

# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = "R/analysis/condition_model_cf_cache/html")

Code for map plots

# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
      axis.text.x = element_text(angle = 90),
      axis.text = element_text(size = 8),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 8),
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
      strip.background = element_rect(fill = "grey95")
      )
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 6),
        strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
        strip.background = element_rect(fill = "gray95"),
        legend.position = c(0.7, 0.02),
        legend.direction = "horizontal"
      )
}

# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2

plot_map_raster <- 
ggplot(swe_coast_proj) + 
  xlim(xmin2, xmax2) +
  ylim(ymin2, ymax2) +
  labs(x = "Longitude", y = "Latitude") +
  geom_sf(size = 0.3) + 
  theme_plot() +
  NULL

plot_map_raster_labels <- 
  plot_map_raster + 
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = 1.9, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = 1.9, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = 1.9, angle = 75)

Read data

# Read the split files and join them
d1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cond_(1_2).csv")
d2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cond_(2_2).csv")

d <- bind_rows(d1, d2)

unique(is.na(d$density_cod))
#> [1] FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE  TRUE

# Calculate standardized variables
d <- d %>% 
  mutate(ln_length_cm = log(length_cm),
         ln_weight_g = log(weight_g),
         year_ct = year - mean(year),
         biomass_her_sc = biomass_her,
         biomass_her_sd_sc = biomass_her_sd,
         biomass_spr_sc = biomass_spr,
         biomass_spr_sd_sc = biomass_spr_sd,
         density_cod_sc = density_cod,
         density_cod_rec_sc = density_cod_rec,
         density_fle_sc = density_fle,
         density_fle_rec_sc = density_fle_rec,
         density_saduria_sc = density_saduria,
         density_saduria_rec_sc = density_saduria_rec,
         depth_sc = depth,
         depth_rec_sc = depth_rec,
         depth_sd_sc = depth_sd,
         oxy_sc = oxy,
         oxy_rec_sc = oxy_rec,
         oxy_sd_sc = oxy_sd,
         temp_sc = temp,
         temp_rec_sc = temp_rec,
         temp_sd_sc = temp_sd,
         year_f = as.factor(year))  %>%
  mutate_at(c("biomass_her_sc", "biomass_her_sd_sc", "biomass_spr_sc", "biomass_spr_sd_sc",
              "density_cod_sc", "density_cod_rec_sc", 
              "density_fle_sc", "density_fle_rec_sc", 
              "density_saduria_sc", "density_saduria_rec_sc", 
              "depth_sc", "depth_rec_sc", "depth_sd_sc",
              "oxy_sc", "oxy_rec_sc", "oxy_sd_sc",
              "temp_sc", "temp_rec_sc", "temp_sd_sc"
              ),
            ~(scale(.) %>% as.vector)) %>% 
  mutate(year = as.integer(year))

unique(is.na(d))
#>      weight_g length_cm  year   lat   lon ices_rect sub_div depth   oxy  temp
#> [1,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [2,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [3,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [4,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [5,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [6,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [7,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#>          X     Y density_cod density_fle density_cod_data density_fle_data
#> [1,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [2,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [3,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [4,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#> [5,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#> [6,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [7,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#>      density_saduria density_saduria_rec density_cod_rec density_fle_rec
#> [1,]           FALSE               FALSE           FALSE           FALSE
#> [2,]           FALSE               FALSE           FALSE           FALSE
#> [3,]           FALSE               FALSE           FALSE           FALSE
#> [4,]           FALSE               FALSE           FALSE           FALSE
#> [5,]           FALSE                TRUE            TRUE            TRUE
#> [6,]           FALSE                TRUE            TRUE            TRUE
#> [7,]           FALSE               FALSE           FALSE           FALSE
#>      depth_rec oxy_rec temp_rec biomass_spr biomass_her density_saduria_sd
#> [1,]     FALSE   FALSE    FALSE       FALSE       FALSE              FALSE
#> [2,]     FALSE   FALSE    FALSE        TRUE        TRUE              FALSE
#> [3,]     FALSE   FALSE    FALSE        TRUE        TRUE              FALSE
#> [4,]     FALSE   FALSE    FALSE       FALSE       FALSE              FALSE
#> [5,]      TRUE    TRUE     TRUE        TRUE        TRUE              FALSE
#> [6,]      TRUE    TRUE     TRUE        TRUE        TRUE              FALSE
#> [7,]     FALSE   FALSE    FALSE        TRUE        TRUE              FALSE
#>      density_cod_sd density_fle_sd depth_sd oxy_sd temp_sd biomass_spr_sd
#> [1,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [2,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [3,]          FALSE          FALSE    FALSE  FALSE   FALSE           TRUE
#> [4,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [5,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [6,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [7,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#>      biomass_her_sd ln_length_cm ln_weight_g year_ct biomass_her_sc
#> [1,]          FALSE        FALSE       FALSE   FALSE          FALSE
#> [2,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [3,]           TRUE        FALSE       FALSE   FALSE           TRUE
#> [4,]          FALSE        FALSE       FALSE   FALSE          FALSE
#> [5,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [6,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [7,]          FALSE        FALSE       FALSE   FALSE           TRUE
#>      biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc density_cod_sc
#> [1,]             FALSE          FALSE             FALSE          FALSE
#> [2,]             FALSE           TRUE             FALSE          FALSE
#> [3,]              TRUE           TRUE              TRUE          FALSE
#> [4,]             FALSE          FALSE             FALSE          FALSE
#> [5,]             FALSE           TRUE             FALSE          FALSE
#> [6,]             FALSE           TRUE             FALSE          FALSE
#> [7,]             FALSE           TRUE             FALSE          FALSE
#>      density_cod_rec_sc density_fle_sc density_fle_rec_sc density_saduria_sc
#> [1,]              FALSE          FALSE              FALSE              FALSE
#> [2,]              FALSE          FALSE              FALSE              FALSE
#> [3,]              FALSE          FALSE              FALSE              FALSE
#> [4,]              FALSE          FALSE              FALSE              FALSE
#> [5,]               TRUE          FALSE               TRUE              FALSE
#> [6,]               TRUE          FALSE               TRUE              FALSE
#> [7,]              FALSE          FALSE              FALSE              FALSE
#>      density_saduria_rec_sc depth_sc depth_rec_sc depth_sd_sc oxy_sc oxy_rec_sc
#> [1,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#> [2,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#> [3,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#> [4,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#> [5,]                   TRUE    FALSE         TRUE       FALSE  FALSE       TRUE
#> [6,]                   TRUE    FALSE         TRUE       FALSE  FALSE       TRUE
#> [7,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#>      oxy_sd_sc temp_sc temp_rec_sc temp_sd_sc year_f
#> [1,]     FALSE   FALSE       FALSE      FALSE  FALSE
#> [2,]     FALSE   FALSE       FALSE      FALSE  FALSE
#> [3,]     FALSE   FALSE       FALSE      FALSE  FALSE
#> [4,]     FALSE   FALSE       FALSE      FALSE  FALSE
#> [5,]     FALSE   FALSE        TRUE      FALSE  FALSE
#> [6,]     FALSE   FALSE        TRUE      FALSE  FALSE
#> [7,]     FALSE   FALSE       FALSE      FALSE  FALSE

unique(is.na(d$density_cod_rec))
#> [1] FALSE  TRUE
unique(is.na(d$density_cod))
#> [1] FALSE

# Drop NA covariates for the correlation plot
d <- d %>% drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
                   density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc, 
                   density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
                   oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc)



# Plot correlation between variables
d_cor <- d %>% dplyr::select("biomass_her_sc", "biomass_her_sd_sc",
                             "biomass_spr_sc", "biomass_spr_sd_sc",
                             "density_cod_sc", "density_cod_rec_sc", 
                             "density_fle_sc", "density_fle_rec_sc", 
                             "density_saduria_sc", "density_saduria_rec_sc", 
                             "depth_sc", "depth_rec_sc",
                             "oxy_sc", "oxy_rec_sc",
                             "temp_sc", "temp_rec_sc")

corr <- round(cor(d_cor), 1)

ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.3))

ggsave("figures/supp/condition/correlation_variables.png", width = 6.5, height = 6.5, dpi = 600)

Read the prediction grids

pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")

pred_grid <- bind_rows(pred_grid1, pred_grid2)

unique(is.na(pred_grid$density_cod))
#> [1] FALSE
unique(is.na(pred_grid$density_cod_rec))
#> [1] FALSE

pred_grid <- pred_grid %>% mutate(year = as.integer(year),
                                  year_f = as.factor(year))

# Scale the variables with respect to data! First calculate means in data
data_means <- d %>%
  dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
                density_cod, density_cod_rec, 
                density_fle, density_fle_rec, 
                density_saduria, density_saduria_rec, 
                depth, depth_rec, depth_sd,
                oxy, oxy_rec, oxy_sd,
                temp, temp_rec, temp_sd) %>%
  mutate_all(~mean(.)) %>%
  distinct(.keep_all = TRUE)

data_stdev <- d %>%
  dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
                density_cod, density_cod_rec, 
                density_fle, density_fle_rec, 
                density_saduria, density_saduria_rec, 
                depth, depth_rec, 
                oxy, oxy_rec, 
                temp, temp_rec
                ) %>%
  mutate_all(~sd(.)) %>%
  distinct(.keep_all = TRUE)

# Before actually scaling, replace the ices-rectangle pelagic values that are NA with the mean across rectangles 
# in the sub division. Replace the sub-division pelagic values that are NA with the mean in the year.
# Note that the sub-division values are not the sum of the rectangles (due to missing rectangles), so 
# I need to calculate a sub-division mean across rectangles within each sub-division
pred_grid <- pred_grid %>% 
  group_by(year) %>% 
  mutate(median_biomass_sprat_across_sub_div = median(biomass_spr_sd, na.rm = TRUE),
         median_biomass_herring_across_sub_div = median(biomass_her_sd, na.rm = TRUE)) %>% 
  ungroup() %>% # Replace sub_divsion NA's with the median across sub_divisions in that year
  mutate(biomass_spr_sd = ifelse(is.na(biomass_spr_sd) == TRUE, median_biomass_sprat_across_sub_div, biomass_spr_sd),
         biomass_her_sd = ifelse(is.na(biomass_her_sd) == TRUE, median_biomass_herring_across_sub_div, biomass_her_sd)) %>% 
  group_by(year, sub_div) %>% 
  mutate(median_biomass_sprat_across_rect_in_sub_div = median(biomass_spr, na.rm = TRUE),
         median_biomass_herring_across_rect_in_sub_div = median(biomass_her, na.rm = TRUE)) %>% 
  ungroup() %>% # Replace rectangle NA's with the median across rectangles in that year and sub-division
  mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, median_biomass_sprat_across_rect_in_sub_div, biomass_spr),
         biomass_her = ifelse(is.na(biomass_her) == TRUE, median_biomass_herring_across_rect_in_sub_div, biomass_her)) %>% 
  # Since I still have some NAs (some sub-divisions do not have a single rectangle in some years), I will will those rectangles 
  # with the sub-division value divided by the number of rectangles in that sub division
  group_by(year, sub_div) %>% 
  mutate(n_rect = length(unique(ices_rect))) %>% 
  ungroup() %>% 
  mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, biomass_spr_sd/n_rect, biomass_spr),
         biomass_her = ifelse(is.na(biomass_her) == TRUE, biomass_her_sd/n_rect, biomass_her))

pred_grid <- pred_grid %>%
  mutate(ln_length_cm = log(1),
         year_ct = year - mean(year),
         biomass_her_sc = (biomass_her - data_means$biomass_her) / data_stdev$biomass_her,
         biomass_her_sd_sc = (biomass_her_sd - data_means$biomass_her_sd) / data_stdev$biomass_her_sd,
         biomass_spr_sc = (biomass_spr - data_means$biomass_spr) / data_stdev$biomass_spr,
         biomass_spr_sd_sc = (biomass_spr_sd - data_means$biomass_spr_sd) / data_stdev$biomass_spr_sd,
         density_cod_sc = (density_cod - data_means$density_cod) / data_stdev$density_cod,
         density_cod_rec_sc = (density_cod_rec - data_means$density_cod_rec) / data_stdev$density_cod_rec,
         density_fle_sc = (density_fle - data_means$density_fle) / data_stdev$density_fle,
         density_fle_rec_sc = (density_fle_rec - data_means$density_fle_rec) / data_stdev$density_fle_rec,
         density_saduria_sc = (density_saduria - data_means$density_saduria) / data_stdev$density_saduria,
         density_saduria_rec_sc = (density_saduria_rec - data_means$density_saduria_rec) / data_stdev$density_saduria_rec,
         depth_sc = (depth - data_means$depth) / data_stdev$depth,
         depth_rec_sc = (depth_rec - data_means$depth_rec) / data_stdev$depth_rec,
         oxy_sc = (oxy - data_means$oxy) / data_stdev$oxy,
         oxy_rec_sc = (oxy_rec - data_means$oxy_rec) / data_stdev$oxy_rec,
         temp_sc = (temp - data_means$temp) / data_stdev$temp,
         temp_rec_sc = (temp_rec - data_means$temp_rec) / data_stdev$temp_rec,
         )

# Plot on map:
ggplot(pred_grid2, aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

pred_grid %>% 
  drop_na() %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

pred_grid %>% 
  drop_na(biomass_spr_sd) %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

unique(is.na(pred_grid))
#>          X     Y depth  year   oxy  temp   lon   lat sub_div ices_rect
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE   FALSE     FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE   FALSE     FALSE
#>      density_saduria biomass_spr biomass_her biomass_spr_sd biomass_her_sd
#> [1,]           FALSE       FALSE       FALSE          FALSE          FALSE
#> [2,]           FALSE       FALSE       FALSE          FALSE          FALSE
#>      density_cod density_fle depth_rec temp_rec oxy_rec density_cod_rec
#> [1,]       FALSE       FALSE     FALSE    FALSE   FALSE           FALSE
#> [2,]       FALSE       FALSE     FALSE    FALSE   FALSE           FALSE
#>      density_fle_rec density_saduria_rec depth_sd temp_sd oxy_sd density_cod_sd
#> [1,]           FALSE               FALSE    FALSE   FALSE  FALSE          FALSE
#> [2,]           FALSE               FALSE    FALSE   FALSE  FALSE          FALSE
#>      density_fle_sd density_saduria_sd year_f
#> [1,]          FALSE              FALSE  FALSE
#> [2,]          FALSE              FALSE  FALSE
#>      median_biomass_sprat_across_sub_div median_biomass_herring_across_sub_div
#> [1,]                               FALSE                                 FALSE
#> [2,]                               FALSE                                 FALSE
#>      median_biomass_sprat_across_rect_in_sub_div
#> [1,]                                       FALSE
#> [2,]                                        TRUE
#>      median_biomass_herring_across_rect_in_sub_div n_rect ln_length_cm year_ct
#> [1,]                                         FALSE  FALSE        FALSE   FALSE
#> [2,]                                          TRUE  FALSE        FALSE   FALSE
#>      biomass_her_sc biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc
#> [1,]          FALSE             FALSE          FALSE             FALSE
#> [2,]          FALSE             FALSE          FALSE             FALSE
#>      density_cod_sc density_cod_rec_sc density_fle_sc density_fle_rec_sc
#> [1,]          FALSE              FALSE          FALSE              FALSE
#> [2,]          FALSE              FALSE          FALSE              FALSE
#>      density_saduria_sc density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc
#> [1,]              FALSE                  FALSE    FALSE        FALSE  FALSE
#> [2,]              FALSE                  FALSE    FALSE        FALSE  FALSE
#>      oxy_rec_sc temp_sc temp_rec_sc
#> [1,]      FALSE   FALSE       FALSE
#> [2,]      FALSE   FALSE       FALSE

pred_grid %>% 
  drop_na(biomass_spr_sc, biomass_spr_sd_sc, biomass_her_sc, biomass_her_sd_sc) %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

Make spde mesh

spde <- make_mesh(d, xy_cols = c("X", "Y"),
                  n_knots = 100,  
                  type = "kmeans", seed = 42)

# Plot and save spde
png(file = "figures/supp/condition/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()

Fit models

Non spatial model to get the average length weight relationship. The main model’s response variable is the ratio of the observed vs this predicted average weight.

ptm <- proc.time()
mnull <- sdmTMB(
  formula = ln_weight_g ~ ln_length_cm,
  data = d,
  time = NULL,
  mesh = spde,
  family = student(link = "identity", df = 5),
  spatiotemporal = "off",
  spatial = "off",
  silent = TRUE,
  reml = FALSE
  )
proc.time() - ptm
#>    user  system elapsed 
#>   1.945   0.263   2.765

# Extract average a and b across the total dataset
a <- tidy(mnull)[1, "estimate"]
b <- tidy(mnull)[2, "estimate"]
a
#> [1] -4.638063
b
#> [1] 2.983866
d$weight_g_avg_pred <- exp(a)*d$length_cm^b

plot(d$weight_g ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$weight_g); abline(a = 0, b = 1, col = "red")

d$le_cren <- d$weight_g / d$weight_g_avg_pred
d$log_le_cren <- log(d$le_cren)
hist(d$le_cren)

Compare breakpoint vs linear oxygen effect in simplified model

mnull_oxy_bp <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + breakpt(oxy_sc) + oxy_rec_sc,
  time_varying = NULL,
  data = d,
  time = "year",
  mesh = spde,
  family = student(link = "identity", df = 5), 
  spatiotemporal = "iid",
  spatial = "on",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1))

mnull_oxy <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + oxy_sc + oxy_rec_sc,
  time_varying = NULL,
  data = d,
  time = "year",
  mesh = spde,
  family = student(link = "identity", df = 5), 
  spatiotemporal = "iid",
  spatial = "on",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1))

AIC(mnull_oxy_bp, mnull_oxy)

min(AIC(mnull_oxy_bp), AIC(mnull_oxy)) 
#> [1] -144518.8

sanity(mnull_oxy_bp)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
sanity(mnull_oxy)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large

tidy(mnull_oxy)

tidy(mnull_oxy_bp)


nd_bp_simple <- data.frame(
  oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100), 
  year = 2019L,
  year_f = as.factor(2019),
  oxy_rec_sc = 0
)

p_bp_simple <- predict(mnull_oxy_bp, newdata = nd_bp_simple,
                       se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))

ggplot(p_bp_simple, aes(oxy_sc, est, 
  ymin = est - 1.96 * est_se, 
  ymax = est + 1.96 * est_se)) +
  geom_line() + geom_ribbon(alpha = 0.4)

ggplot(d, aes(oxy_sc, oxy)) + 
  geom_point() + 
  geom_hline(yintercept = (1.0700*sd(d$oxy)) + mean(d$oxy)) +
  geom_vline(xintercept = 1.0700)

Fit main model

ptm <- proc.time()
mfull <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
    biomass_spr_sc + biomass_spr_sd_sc +
    density_cod_sc + density_cod_rec_sc + 
    density_fle_sc + density_fle_rec_sc + 
    density_saduria_sc + density_saduria_rec_sc + 
    depth_sc + depth_rec_sc + 
    oxy_sc + oxy_rec_sc + 
    temp_sc + temp_rec_sc,
    data = d,
  time = "year",
  mesh = spde, 
  family = student(link = "identity", df = 5),
  spatiotemporal = "iid",
  spatial = "on",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1)
  )
proc.time() - ptm
#>    user  system elapsed 
#> 721.928  20.345 747.627

tidy(mfull, conf.int = TRUE) %>% arrange(desc(abs(estimate)))

tidy(mfull, effects = "ran_par", conf.int = TRUE)

# Check model convergence
sanity(mfull)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large

mcmc_res <- residuals(mfull, type = "mle-mcmc", mcmc_iter = 5000, mcmc_warmup = 2000)

png(file = "figures/supp/condition/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()

d$residuals <- mcmc_res[, 1]

Calculate approximate r^2 manually using a modifed r2.sdmTMB

b_fe <- tidy(mfull)
b_fe <- stats::setNames(b_fe$estimate, b_fe$term)
X <- mfull$tmb_data$X_ij
X <- matrix(unlist(X), ncol = length(b_fe))

VarF <- var(as.vector(b_fe %*% t(X))) # variance from fixed-effects
 
b_ran <- tidy(mfull, "ran_par")
sigma_O <- b_ran$estimate[b_ran$term == "sigma_O"] # spatial standard deviation
sigma_E <- b_ran$estimate[b_ran$term == "sigma_E"] # spatiotemporal standard deviation
VarO <- sigma_O^2 # spatial variance
VarE <- sigma_E^2 # spatiotemporal variance
 
# Calculate obs - pred
d_r2 <- d
d_r2$pred <- predict(mfull)
d_r2$resid_manual <- d_r2$le_cren - d_r2$pred$est
#hist(d_r2$resid_manual)

VarR <- var(as.vector(d_r2$resid_manual)) # "residual" variance

# Marginal R2s
denominator <- VarF + VarO + VarE + VarR

marg <- VarF/denominator

out <- data.frame(
    marginal = marg,
    partial_spatial = VarO/denominator,
    partial_spatiotemporal = VarE/denominator,
    conditional = (VarF + VarO + VarE)/denominator,
    all_random = (VarO + VarE)/denominator,
    marginal_random_ratio = marg / ((VarO + VarE)/denominator),
    random_marginal_ratio = ((VarO + VarE)/denominator) / marg
  )

out


# What if I skip the year effect?
b_fe2 <- tidy(mfull) %>% filter(!grepl('year', term))
#> filter: removed 27 rows (63%), 16 rows remaining
b_fe2 <- stats::setNames(b_fe2$estimate, b_fe2$term)

X2 <- mfull$tmb_data$X_ij
X2 <- matrix(unlist(X2), ncol = length(b_fe))[, 28:43] # skip year columns
VarF2 <- var(as.vector(b_fe2 %*% t(X2))) # variance from fixed-effects

denominator3 <- VarF2 + VarO + VarE + VarR

marg3 <- VarF2/denominator3
marg3
#> [1] 0.05794511

Plot residuals

# Residuals on map
plot_map_raster +
  geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals)) +
  scale_colour_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  labs(color = "residuals") +
  NULL


ggsave("figures/supp/condition/spatial_resid.png", width = 6.5, height = 6.5, dpi = 600)

# Residuals vs length and year
ggplot(d, aes(x = length_cm, y = residuals, color = residuals)) +
  geom_point(size = 0.1, alpha = 1) +
  geom_hline(yintercept = 0, size = 0.5, color = "black", linetype = 2, alpha = 0.5) +
  labs(x = "length [cm]", y = "Residuals", color = "") +
  stat_smooth(color = "black", size = 0.5) +
  facet_wrap(~year, ncol = 5) +
  scale_color_gradient2() +
  theme_facet_map()


ggsave("figures/supp/condition/residuals_vs_length_and_year.png", width = 6.5, height = 6.5, dpi = 600)

ggplot(d, aes(year, length_cm, z = residuals)) +
  stat_summary_2d(bins = 40) +
  scale_fill_gradient2()


# ggsave("figures/supp/condition/residuals_vs_length_and_year2.png", width = 6.5, height = 6.5, dpi = 600)

Annual index using `get_index_sims` method to avoid slow bias correction

ncells <- filter(pred_grid, year == max(pred_grid$year)) %>% nrow()
#> filter: removed 240,214 rows (96%), 9,239 rows remaining

pred_avg_sim <- predict(mfull, newdata = pred_grid, nsim = 100)

index_avg_sim <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_sims <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)),
                                 return_sims = TRUE) # This is for calculation on samples!
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

write.csv(index_avg_sim, "output/index_avg_sim.csv")

Evaluate sensitivity to using a pred grid with a subset excluding low density areas

# Now do the same but excluding the areas not sampled in the pred grid
ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 20 & depth < 120) %>% nrow()
#> filter: removed 242,779 rows (97%), 6,674 rows remaining
pred_grid_sub <- filter(pred_grid, depth > 20 & depth < 120)
#> filter: removed 69,255 rows (28%), 180,198 rows remaining

pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 100)

index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(1/ncells_sub, nrow(pred_avg_sim_sub)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

# Combine and plot & compare
index_avg_sim_comp <- bind_rows(mutate(index_avg_sim, area = "full"),
                                mutate(index_avg_sim_sub, area = "subset"))
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'area' (character) with one unique value and 0% NA

ggplot(index_avg_sim_comp, aes(year, log(est), ymin = log(lwr), ymax = log(upr), color = area)) +
  geom_line() +
  geom_ribbon(alpha = 0.4)

ggplot(index_avg_sim_comp, aes(year, est, ymin = lwr, ymax = upr, color = area)) +
  geom_line() +
  geom_ribbon(alpha = 0.4)

Difference in Le Cren condition factor between the end and start of time period using sims

colnames(pred_grid)
#>  [1] "X"                                            
#>  [2] "Y"                                            
#>  [3] "depth"                                        
#>  [4] "year"                                         
#>  [5] "oxy"                                          
#>  [6] "temp"                                         
#>  [7] "lon"                                          
#>  [8] "lat"                                          
#>  [9] "sub_div"                                      
#> [10] "ices_rect"                                    
#> [11] "density_saduria"                              
#> [12] "biomass_spr"                                  
#> [13] "biomass_her"                                  
#> [14] "biomass_spr_sd"                               
#> [15] "biomass_her_sd"                               
#> [16] "density_cod"                                  
#> [17] "density_fle"                                  
#> [18] "depth_rec"                                    
#> [19] "temp_rec"                                     
#> [20] "oxy_rec"                                      
#> [21] "density_cod_rec"                              
#> [22] "density_fle_rec"                              
#> [23] "density_saduria_rec"                          
#> [24] "depth_sd"                                     
#> [25] "temp_sd"                                      
#> [26] "oxy_sd"                                       
#> [27] "density_cod_sd"                               
#> [28] "density_fle_sd"                               
#> [29] "density_saduria_sd"                           
#> [30] "year_f"                                       
#> [31] "median_biomass_sprat_across_sub_div"          
#> [32] "median_biomass_herring_across_sub_div"        
#> [33] "median_biomass_sprat_across_rect_in_sub_div"  
#> [34] "median_biomass_herring_across_rect_in_sub_div"
#> [35] "n_rect"                                       
#> [36] "ln_length_cm"                                 
#> [37] "year_ct"                                      
#> [38] "biomass_her_sc"                               
#> [39] "biomass_her_sd_sc"                            
#> [40] "biomass_spr_sc"                               
#> [41] "biomass_spr_sd_sc"                            
#> [42] "density_cod_sc"                               
#> [43] "density_cod_rec_sc"                           
#> [44] "density_fle_sc"                               
#> [45] "density_fle_rec_sc"                           
#> [46] "density_saduria_sc"                           
#> [47] "density_saduria_rec_sc"                       
#> [48] "depth_sc"                                     
#> [49] "depth_rec_sc"                                 
#> [50] "oxy_sc"                                       
#> [51] "oxy_rec_sc"                                   
#> [52] "temp_sc"                                      
#> [53] "temp_rec_sc"

# Plot the sims
ggplot(index_avg_sims, aes(y = as.factor(year), x = .value)) + 
  geom_density_ridges_gradient(rel_min_height = 0.01, alpha = 0.8, color = "grey10", size = 0.2,
                               quantile_lines = TRUE, quantiles = 0.5) + 
  coord_flip()
#> Picking joint bandwidth of 0.00302


# index_avg_sims %>% 
#   filter(year %in% c(1993, 2019)) %>% 
#   ggplot(., aes(as.factor(year), .value)) +
#   geom_violin(width = 1.3)

# Now reshape data again to calculate row-wise differences
year_diff <- index_avg_sims %>% 
  filter(year %in% c(1993, 2019)) %>% 
  #mutate(.value = log(.value)) %>% 
  pivot_wider(names_from = year, values_from = .value) %>% 
  mutate("2019-1993" = `2019` - `1993`) %>% 
  pivot_longer(2:4) %>% 
  ungroup()
#> filter: removed 2,500 rows (93%), 200 rows remaining
#> pivot_wider: reorganized (year, .value) into (1993, 2019) [was 200x3, now 100x3]
#> mutate: new variable '2019-1993' (double) with 100 unique values and 0% NA
#> pivot_longer: reorganized (1993, 2019, 2019-1993) into (name, value) [was 100x4, now 300x3]
#> ungroup: no grouping variables

# Plot these
p1 <- ggplot(year_diff) +
  geom_density(data = filter(year_diff, !name == "2019-1993"),
               aes(value, fill = name), alpha = 0.6) +
  scale_fill_brewer(palette = "Dark2") +
  coord_cartesian(expand = 0)
#> filter: removed 100 rows (33%), 200 rows remaining

p2 <- ggplot(year_diff) +
  geom_density(data = filter(year_diff, name == "2019-1993"),
               aes(value), alpha = 0.6, fill = "grey") +
  geom_vline(xintercept = 0, linetype = 2) 
#> filter: removed 200 rows (67%), 100 rows remaining

p1/p2


# Calculate quantiles for each level (to go in the results section)
year_diff %>% 
  group_by(name) %>% 
  summarise(median = quantile(value, probs = 0.5),
            upr97.5 = quantile(value, probs = 0.975),
            lwr2.5 = quantile(value, probs = 0.025))
#> group_by: one grouping variable (name)
#> summarise: now 3 rows and 4 columns, ungrouped


# Percent change in condition factor
year_diff %>% 
  group_by(name) %>% 
  pivot_wider(names_from = name, values_from = value) %>%
  mutate(perc_change = ((`2019`-`1993`)/`1993`)*100) %>% 
  ungroup() %>% 
  summarise(lwr2.5 = quantile(perc_change, probs = 0.025),
            median = quantile(perc_change, probs = 0.5),
            upr97.5 = quantile(perc_change, probs = 0.975))
#> group_by: one grouping variable (name)
#> pivot_wider: reorganized (name, value) into (1993, 2019, 2019-1993) [was 300x3, now 100x4]
#> mutate: new variable 'perc_change' (double) with 100 unique values and 0% NA
#> ungroup: no grouping variables
#> summarise: now one row and 3 columns, ungrouped

Simulated annual condition factor by sub-division

# From these models, predict annual condition factor
# I will use the prediction grid with ALL covariates set to 0, incl. depth, temp and oxygen
# Grabbing the number of cells to help with calculating the average
ncells24 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 24) %>% nrow()
ncells25 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 25) %>% nrow()
ncells26 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 26) %>% nrow()
ncells27 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 27) %>% nrow()
ncells28 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 28) %>% nrow()

# Use the `area` argument here to turn the total into an average by giving it one over the number of cells
pred_avg_24_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 24), sims = 100)
#> Warning: The `sims` argument of `predict.sdmTMB()` is deprecated as of sdmTMB 0.0.21.
#> Please use the `nsim` argument instead.
pred_avg_25_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 25), sims = 100)
pred_avg_26_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 26), sims = 100)
pred_avg_27_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 27), sims = 100)
pred_avg_28_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 28), sims = 100)

index_avg_24_sim <- get_index_sims(pred_avg_24_sim, area = rep(1/ncells24, nrow(pred_avg_24_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_25_sim <- get_index_sims(pred_avg_25_sim, area = rep(1/ncells25, nrow(pred_avg_25_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_26_sim <- get_index_sims(pred_avg_26_sim, area = rep(1/ncells26, nrow(pred_avg_26_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_27_sim <- get_index_sims(pred_avg_27_sim, area = rep(1/ncells27, nrow(pred_avg_27_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_28_sim <- get_index_sims(pred_avg_28_sim, area = rep(1/ncells28, nrow(pred_avg_28_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

# Combine and plot
index_avg_24_sim <- index_avg_24_sim %>% mutate(sub_div = 24)
index_avg_25_sim <- index_avg_25_sim %>% mutate(sub_div = 25)
index_avg_26_sim <- index_avg_26_sim %>% mutate(sub_div = 26)
index_avg_27_sim <- index_avg_27_sim %>% mutate(sub_div = 27)
index_avg_28_sim <- index_avg_28_sim %>% mutate(sub_div = 28)

div_index <- bind_rows(index_avg_24_sim, index_avg_25_sim, index_avg_26_sim, index_avg_27_sim, index_avg_28_sim)

Condition factor over time and individual distributions

div_index2 <- bind_rows(mutate(div_index, sub_div = as.factor(sub_div)),
                        mutate(index_avg_sim, sub_div = as.factor("Total")))

pal <- c(brewer.pal(n = 5, name = "Dark2"), "grey10")

# Remember to log now that we use sims
cond_index <- ggplot(div_index2, aes(year, est, color = sub_div, linetype = sub_div)) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  labs(y = "Predicted Le Cren's condition factor", x = "Year") +
  geom_line(alpha = 0.8) +
  geom_point(alpha = 0) + 
  scale_color_manual(values = pal, name = "ICES subdivision") +
  scale_linetype_manual(values = c(1,1,1,1,1,0)) +
  geom_point(data = index_avg_sim, aes(year, est), inherit.aes = FALSE, 
             size = 2.5, shape = 21, fill = "gray1", color = "white") +
  geom_errorbar(data = index_avg_sim, inherit.aes = FALSE, 
                aes(x = year, ymax = upr, ymin = lwr),
                width = 0, alpha = 0.8, color = "gray1", size = 0.7) +
  theme(axis.title = element_text(size = 9),
        legend.position = "bottom",
        legend.background = element_blank(),
        legend.title = element_text(size = 9),
        legend.spacing.x = unit(0, 'cm'),
        plot.margin = unit(c(0.4, 0.4, 0.4, 0), "cm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  guides(linetype = "none",
         color = guide_legend(title.position = "top", title.hjust = 0.5, nrow = 1,
                              keywidth = 0.1, keyheight = 0.1, default.unit = "inch",
                              override.aes = list(linetype = c(1,1,1,1,1,0),
                                                  shape = c(NA,NA,NA,NA,NA,16),
                                                  alpha = rep(0.8, 6),
                                                  color = pal))) + 
  NULL

cond_index


# Ridge-plot of data
data_ridge <- ggplot(d, aes(x = le_cren, y = fct_rev(as.factor(year)), fill = factor(stat(quantile)))) +
  stat_density_ridges(
    rel_min_height = .01,
    scale = 3,
    geom = "density_ridges_gradient",
    calc_ecdf = TRUE,
    alpha = 0.1,
    quantiles = c(0.25, 0.75)) +
  scale_fill_manual(
    name = "Probability",
    values = c("grey50", "gray70", "grey90"),
    labels = c("[0,0.25]", "[0.25,0.75]", "[0.75,1]")) +
  xlim(0.65, 1.55) + 
  labs(x = "Le Cren's condidition factor",
       y = "Year") + 
  guides(fill = guide_legend(nrow = 1, title.position = "top", title.hjust = 0.5)) + 
  theme(axis.title = element_text(size = 9),
        legend.position = "bottom",
        legend.key.size = unit(0.5, "line"),
        legend.background = element_blank(),
        legend.title = element_text(size = 9),
        legend.spacing.x = unit(0.1, 'cm'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        plot.margin = unit(c(0.4, 0.8, 0.4, 0), "cm")) +
  NULL

data_ridge + cond_index + plot_annotation(tag_levels = 'A') +
   plot_layout(widths = c(1, 2.2))
#> Warning: Removed 926 rows containing non-finite values (stat_density_ridges).


ggsave("figures/cond_index.pdf", width = 17, height = 17, units = "cm")
#> Warning: Removed 926 rows containing non-finite values (stat_density_ridges).


## Calculate delta_fulton and relate to delta_oxygen, by subdivision
div_index_delta_cond <- div_index2 %>% 
  group_by(sub_div) %>% 
  filter(year %in% c(1993, 1994, 1995, 2017, 2018, 2019)) %>% 
  dplyr::select(year, log_est, sub_div) %>% 
  pivot_wider(names_from = year, values_from = log_est) %>% 
  mutate(start_cond = mean(`1993`, `1994`, `1995`),
         end_cond = mean(`2017`, `2018`, `2019`),
         delta_cond = end_cond - start_cond) %>% 
  ungroup() %>% 
  dplyr::select(sub_div, delta_cond)

write.csv(div_index_delta_cond, "output/div_index_delta_cond.csv")

Effect sizes

# What is one standard deviation of oxygen?
sd(d$oxy)
#> [1] 1.846734

# Extract random and fixed coefficients from the full model
mfull_est <- bind_rows(tidy(mfull, effects = "ran_par", conf.int = TRUE) %>%
                         filter(term %in% c("sigma_O", "sigma_E")),
                       
                       tidy(mfull, effects = "fixed", conf.int = TRUE)  %>% 
                         filter(!grepl('year', term))) %>%
  
  mutate(term = factor(term))

mfull_est <- mfull_est %>% 
  mutate(group = "Herring",
         group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
         group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
         group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
         group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
         group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
         group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
         group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
         group = ifelse(term == "sigma_E", "Random", group),
         group = ifelse(term == "sigma_O", "Random", group)
         )

mfull_est <- mfull_est %>% 
  mutate(scale = "Sub-division",
         scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
                                    "density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
         scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
                                    "depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
         scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
         scale = ifelse(term == "sigma_O", "Spatial/spatiotemporal s.d.", scale)
         )

# Plot effects
# This is the order changed labels should go in!
levels(mfull_est$term)
#>  [1] "biomass_her_sc"         "biomass_her_sd_sc"      "biomass_spr_sc"        
#>  [4] "biomass_spr_sd_sc"      "density_cod_rec_sc"     "density_cod_sc"        
#>  [7] "density_fle_rec_sc"     "density_fle_sc"         "density_saduria_rec_sc"
#> [10] "density_saduria_sc"     "depth_rec_sc"           "depth_sc"              
#> [13] "oxy_rec_sc"             "oxy_sc"                 "sigma_E"               
#> [16] "sigma_O"                "temp_rec_sc"            "temp_sc"

# I want the random effects to be dark gray...
pal <- brewer.pal(n = 8, name = "Dark2")
pal2 <- c(pal[1:5], "black", pal[6:8])

# Sort the terms so that random effects are at the top...
mfull_est <- mfull_est %>% 
  mutate(term2 = ifelse(term %in% c("sigma_E", "sigma_O"), 2, 1))

ggplot(mfull_est, aes(reorder(term, term2), estimate, color = group, fill = group, shape = scale)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray40", alpha = 0.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
  geom_point(size = 2.5) +
  scale_color_manual(values = pal2, name = "") +
  scale_shape_manual(values = c(15, 19, 21, 17), name = "") +
  scale_fill_manual(values = rep("white", 9), name = "") +
  labs(x = "", y = "Standardized coefficient") +
  scale_x_discrete(breaks = levels(mfull_est$term),
                   labels = c(expression(Herring[rec]),
                              expression(Herring[sub]),
                              expression(Sprat[rec]),
                              expression(Sprat[sub]),
                              expression(Cod[rec]),
                              expression(Cod[haul]),
                              expression(Flounder[rec]),
                              expression(Flounder[haul]),
                              expression(Saduria~entomon[rec]),
                              expression(Saduria~entomon[haul]),
                              expression(Depth[rec]),
                              expression(Depth[haul]),
                              expression(Oxygen[rec]),
                              expression(Oxygen[haul]),
                              expression(sigma[E]),
                              expression(sigma[O]),
                              expression(Temp[rec]),
                              expression(Temp[haul])
                              )) +
  coord_flip() +
  theme(plot.margin = unit(c(0.4, 0.4, 0.4, 0), "cm"),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  guides(color = "none", fill = "none", shape = guide_legend(ncol = 2)) + 
  NULL

  
ggsave("figures/effect_size.pdf", width = 17, height = 17, units = "cm")


# Just a test to see the labels were alright
#ggplot(mfull_est, aes(reorder(term, estimate), estimate)) +
# ggplot(mfull_est, aes(term, estimate)) +
#   geom_hline(yintercept = 0, linetype = 2, color = "red", alpha = 0.5) +
#   geom_point(size = 2) +
#   labs(x = "", y = "Standardized coefficient") +
#   coord_flip()

# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est, "output/mfull_est.csv")

Plot in space

pred_map <- predict(mfull, newdata = pred_grid)

pal <- wes_palette("Zissou1", 21, type = "continuous")

plot_map_raster +
  geom_raster(data = filter(pred_map, depth < 120 & depth > 10), aes(x = X*1000, y = Y*1000, fill = exp(est))) +
  scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
  facet_wrap(~ year, ncol = 5) +
  theme_facet_map() +
  NULL


ggsave("figures/supp/condition/all_years_condition_map_covars.png", width = 6.5, height = 6.5, dpi = 600)

# And a smaller map for selected years
plot_map_raster_labels +
  geom_raster(data = filter(pred_map, year %in% c(1994, 2001, 2008, 2018) & depth < 120 & depth > 10),
              aes(x = X*1000, y = Y*1000, fill = exp(est))) +
  scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") + 
  facet_wrap(~ year, ncol = 2) +
  theme_plot() +
  NULL


ggsave("figures/condition_map_covars.pdf", width = 17, height = 17, units = "cm")



# Plot correlation between random effects and covariates

# Drop NA covariates for the correlation plot
# colnames(pred_map)
# 
# pred_map_corrplot <- pred_map %>%
#   drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
#           density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc, 
#           density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
#           oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc, omega_s, epsilon_st) %>% 
#   dplyr::select(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
#           density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc, 
#           density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
#           oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc, omega_s, epsilon_st)
# 
# corr_re <- round(cor(pred_map_corrplot), 2)
# 
# ggcorrplot(corr_re, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
#   theme_classic() + 
#   theme(axis.text.x = element_text(angle = 90, vjust = 0.3))
# 
# ggsave("figures/supp/condition/correlation_variables_random_effects.png", width = 6.5, height = 6.5, dpi = 600)

Plot spatial and spatiotemporal random effects:

# Plot spatiotemporal random effect
plot_map_raster +
  geom_raster(data = pred_map, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  ggtitle("Spatiotemporal random effects") +
  NULL


ggsave("figures/supp/condition/epsilon_st_map.png", width = 6.5, height = 6.5, dpi = 600)

# Plot spatial random effect
plot_map_raster +
  geom_raster(data = filter(pred_map, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  theme_plot() + 
  ggtitle("Spatial random effects") +
  NULL


ggsave("figures/supp/condition/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)

Calculate the “spatial trend” from the estimates:

# Fit a linear model to each prediction grid of the estimate over time
# https://community.rstudio.com/t/extract-slopes-by-group-broom-dplyr/2751/7
time_slopes_by_year <- pred_map %>%
  mutate(id = paste(X, Y, sep = "_")) %>%
  split(.$id) %>%
  purrr::map(~lm(est ~ year, data = .x)) %>%
  purrr::map_df(broom::tidy, .id = 'id') %>%
  filter(term == 'year')

# Plot the slopes
time_slopes_by_year2 <- time_slopes_by_year %>%
  separate(id, c("X", "Y"), sep = "_") %>%
  mutate(X = as.numeric(X),
         Y = as.numeric(Y))

plot_map_raster +
  geom_raster(data = time_slopes_by_year2, aes(x = X * 1000, y = Y * 1000, fill = estimate)) +
  scale_fill_gradient2(midpoint = 0) +
  facet_wrap(~ term, ncol = 5) +
  theme_plot() +
  ggtitle("Time slopes by each pred grid") +
  NULL


ggsave("figures/supp/condition/spatial_trend_condition.png", width = 6.5, height = 6.5, dpi = 600)

Visualize marginal effects

Depth

# Prepare prediction data frame
nd_dep <- data.frame(depth_sc = seq(min(d$depth_sc), max(d$depth_sc), length.out = 100))

nd_dep <- nd_dep %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_dep <- predict(mfull, newdata = nd_dep, se_fit = TRUE, re_form = NA)

mar_dep <- ggplot(p_margin_dep, aes(depth_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Depth[haul]))

Oxygen

# Prepare prediction data frame
nd_oxy <- data.frame(oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100))

nd_oxy <- nd_oxy %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         #oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_oxy <- predict(mfull, newdata = nd_oxy, se_fit = TRUE, re_form = NA)

mar_oxy <- ggplot(p_margin_oxy, aes(oxy_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se)))+
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Oxygen[haul]))

Temperature

# Prepare prediction data frame
nd_tem <- data.frame(temp_sc = seq(min(d$temp_sc), max(d$temp_sc), length.out = 100))

nd_tem <- nd_tem %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         #temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_tem <- predict(mfull, newdata = nd_tem, se_fit = TRUE, re_form = NA)

mar_temp <- ggplot(p_margin_tem, aes(temp_sc, exp(est),
                                     ymin = exp(est - 1.96 * est_se),
                                     ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Temperature[haul]))

Cod

# Prepare prediction data frame
nd_cod <- data.frame(density_cod_sc = seq(min(d$density_cod_sc), max(d$density_cod_sc), length.out = 100))

nd_cod <- nd_cod %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         #density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_cod <- predict(mfull, newdata = nd_cod, se_fit = TRUE, re_form = NA)

mar_cod <- ggplot(p_margin_cod, aes(density_cod_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Cod[haul]))

Sprat

# Prepare prediction data frame
nd_spr <- data.frame(biomass_spr_sd_sc = seq(min(d$biomass_spr_sd_sc), max(d$biomass_spr_sd_sc), length.out = 100))

nd_spr <- nd_spr %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         #biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_spr <- predict(mfull, newdata = nd_spr, se_fit = TRUE, re_form = NA)

mar_spr <- ggplot(p_margin_spr, aes(biomass_spr_sd_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  #coord_cartesian(ylim = c(-4.65, -4.4)) + 
  theme(aspect.ratio = 1) + 
  labs(x = expression(Sprat[sub]))

Saduria

# Prepare prediction data frame
nd_sad <- data.frame(density_saduria_sc = seq(min(d$density_saduria_sc), max(d$density_saduria_sc), length.out = 100))

nd_sad <- nd_sad %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         #density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_sad <- predict(mfull, newdata = nd_sad, se_fit = TRUE, re_form = NA)

mar_sad <- ggplot(p_margin_sad, aes(density_saduria_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  #coord_cartesian(ylim = c(-4.75, -4.4)) + 
  theme(aspect.ratio = 1) + 
  labs(x = expression(Saduria~entomon[haul]))

Plot together

# mar_dep + mar_oxy + mar_temp + mar_cod + mar_spr + mar_sad
#   plot_annotation(tag_levels = 'A') 
# 
# ggsave("figures/supp/condition/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)

p_margin_dep2 <- p_margin_dep %>%
  mutate(variable = "Depth (haul)") %>% 
  dplyr::select(est, est_se, depth_sc, variable) %>% 
  rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_oxy2 <- p_margin_oxy %>%
  mutate(variable = "Oxygen (haul)") %>% 
  dplyr::select(est, est_se, oxy_sc, variable) %>% 
  rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_tem2 <- p_margin_tem %>%
  mutate(variable = "Temperature (haul)") %>% 
  dplyr::select(est, est_se, temp_sc, variable) %>% 
  rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_cod2 <- p_margin_cod %>%
  mutate(variable = "Cod (haul)") %>% 
  dplyr::select(est, est_se, density_cod_sc, variable) %>% 
  rename(value = density_cod_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_spr2 <- p_margin_spr %>%
  mutate(variable = "Sprat (sub)") %>% 
  dplyr::select(est, est_se, biomass_spr_sd_sc, variable) %>% 
  rename(value = biomass_spr_sd_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_sad2 <- p_margin_sad %>%
  mutate(variable = "S. entomon (haul)") %>%
  dplyr::select(est, est_se, density_saduria_sc, variable) %>% 
  rename(value = density_saduria_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

margs <- bind_rows(p_margin_dep2, p_margin_oxy2, p_margin_tem2, 
                   p_margin_cod2, p_margin_spr2, p_margin_sad2)

ggplot(margs, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
                  ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  facet_wrap(~variable, scales = "free_x") +
  labs(y = "Le Cren's condidition factor",
       x = "Scaled value") +
  theme_plot() +
  theme(axis.text.x = element_text(angle = 0)) +
  NULL


ggsave("figures/supp/condition/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)

Condition model

Max Lindmark, Sean C. Andersson, Mayya Gogina and Michele Casini

21 October, 2022

Fit condition model with environmental and biological predictors

Load packages

Code for map plots

Read data

Read the prediction grids

Make spde mesh

Fit models

Non spatial model to get the average length weight relationship. The main model’s response variable is the ratio of the observed vs this predicted average weight.

Compare breakpoint vs linear oxygen effect in simplified model

Fit main model

Calculate approximate r^2 manually using a modifed r2.sdmTMB

Plot residuals

Annual index using `get_index_sims` method to avoid slow bias correction

Difference in Le Cren condition factor between the end and start of time period using sims

Simulated annual condition factor by sub-division

Condition factor over time and individual distributions

Effect sizes

Plot in space

Plot spatial and spatiotemporal random effects:

Calculate the “spatial trend” from the estimates:

Visualize marginal effects

Depth

Oxygen

Temperature

Cod

Sprat

Saduria

Plot together

Condition model

Max Lindmark, Sean C. Andersson, Mayya Gogina and Michele Casini

21 October, 2022

Fit condition model with environmental and biological predictors

Load packages

Code for map plots

Read data

Read the prediction grids

Make spde mesh

Fit models

Non spatial model to get the average length weight relationship. The main model’s response variable is the ratio of the observed vs this predicted average weight.

Compare breakpoint vs linear oxygen effect in simplified model

Fit main model

Calculate approximate r^2 manually using a modifed r2.sdmTMB

Plot residuals

Annual index using get_index_sims method to avoid slow bias correction

Difference in Le Cren condition factor between the end and start of time period using sims

Simulated annual condition factor by sub-division

Condition factor over time and individual distributions

Effect sizes

Plot in space

Plot spatial and spatiotemporal random effects:

Calculate the “spatial trend” from the estimates:

Visualize marginal effects

Depth

Oxygen

Temperature

Cod

Sprat

Saduria

Plot together

Annual index using `get_index_sims` method to avoid slow bias correction